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CARSON, WW. ABSTRACT 

An investigation was na @RAG uckling of rectangular plates 
with free, simply supported, or clamped edges and loaded in uniform 
compression or pure shear. The finite element method with a rectangu- 
lar, sixteen-degree-of-freedom element was used. Results show that, 
for the same order of error, this element allows much coarser plate 
divisions than a widely used twelve-degree-of-freedom element. An 
order of magnitude reduction in the required order of the assembled 


stiffness and stability coefficient matrices results. 
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LIST OF SYMBOLS 
plate length in x direction; polynomial coefficient 
vector of displacement polynomial coefficients 
transformation matrix 
plate length in y direction 
stability coefficient matrix 
plate flexural rigidity, Eh?/(12(1 - y7 5) 
modulus of elasticity 
vector of external nodal forces and moments 


matrix of functions of x and y, relating vectors 
and z 


|m 


plate thickness 


matrix of functions of x and y, relating vectors 
and r 


| 


identity matrix 
buckling coefficient 


stiffness matrix 


w~ 
il 
7) 
w~ 
+ 


modified stiffness matrix from defining relation K 
vector of functions of bending moments 

bending moment per unit length 

load matrix 


matrix of constants giving relative magnitudes of types 
of loads 


plate element length in x direction 

matrix of functions of Poisson's ratio 

plate element width in y direction 

vector of first derivatives of transverse displacement 
strain energy of bending 


total potential energy 


|O, 


transverse displacement 
work of in-plane forces 
rectangular coordinates 


vector of second derivatives of transverse 
displacement 


clamped edge 

free edge 

simply supported edge 

vector of nodal displacements and slopes 
scalar factor of matrix N* 

eigenvalue 

Poisson's ratio 

transpose of a matrix 

inverse of a matrix 


in the x direction; for moments, about a line 
parallel to the x axis 


partial derivative with respect to x 


pertaining to a plate element 


1. INTRODUCTION 
BACKGROUND 

The Finite Element Method for Stability Problems 

The finite element method provides a numerical means of analyzing 
boundary value field problems. When applied to structural mechanics, 
the method provides a means of deriving a stiffness matrix for an as- 
sembly of structural elements. For stability problems, a stability co- 
efficient matrix is also required. This matrix is used in conjunction 
with the stiffness matrix. 
Development of the Method 

Turner, Clough, Martin, and Topp introduced the finite element 
method in 1956." In 1963, Melosh presented a basis for derivation of 
the stiffness matrix for a thin plate element, using three degrees of 
freedom per Agee” Also in 1963, Gallagher and Padliog first used a 
stability coefficient matrix, for column Pehiges In 1966, Kapur and 
Hartz gave the results of the use of the stability coefficient matrix 
approach for plate buckling pee ie: They used the Melosh plate ele- 
ment. to find the stiffness matrix. While Kapur and Hartz were making 
their study of plate buckling, in 1965 Bogner, Fox, and Schmit presented 
a new stiffness matrix for a plate enenede.° Four degrees of freedom 
per node were used for the new element. 
Significance of Degrees of Freedom 

For the Melosh method of deriving the stiffness matrix of a plate 
element, the three degrees of freedom per node are lateral deflection 
and slope in each of two perpendicular, in-plane directions (w; Wes 
wi This assures compatibility of Lateral deflection of adjacent edges 
of elements. However, complete continuity of slope is not achieved. The 


o 


method of Bogner, Fox, and Schmit includes an additional degree of 
freedom — which does provide continuity of slope. For a rectangu- 
lar plate element, an element stiffness matrix of order sixteen is re- 
quired to represent four degrees of freedom at each of four nodes. For 
three degrees of freedom per node, the order is twelve. Hence, the six- 
teen-degree-of-freedom element stiffness matrix is larger but should 


produce more accurate results because of better continuity of slope.* 


THE PROBLEM INVESTIGATED 

Ob jective 

The objective of the investigation was to determine the consequences 
of the use of the newer, sixteen-degree-of-freedom element stiffness 
matrix in plate buckling problems. Secondarily, some previously unsolved 
cases of plate buckling were investigated to show the versatility of the 
finite element method. 
Limitations 

The investigation was limited to finding buckling coefficients and 
buckling mode vectors for thin, homogeneous, isotropic, rectangular plates 
with free, simply supported, or clamped edge conditions, and with uniform 


compression or pure shear loads. Plates were subdivided into identical 


rectangular elements. 


JUSTIFICATION OF THE INVESTIGATION 


Advantages of the Finite Element Method 


*Other rectangular plate elements with twelve or sixteen degrees of 
freedom may be formed by using different displacement functions. Through- 
out, plate elements with twelve or sixteen degrees of freedom will mean 
the elements discussed above. 
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Plate buckling problems of even moderate complexity require numeri- 
cal procedures. Of the numerical procedures applicable, the finite ele- 
ment method has some unique advantages, as follows: 


ie The capability of analyzing plates of arbitrary shape, such 
as plates with holes. 


2 The capability of analyzing plates made of materials with 
orthotropic properties. 


Sg The capability of introducing edge conditions conveniently, 
especially the free edge. 


4. The capability of synthesizing complex structures from plates 
and other members, such as beams. 


Advantages of Using Sixteen-Degree-of-Freedom Element 





Assembled stiffness and stability coefficient matrices may be very 
large, especially in structural problems which require synthesis of 
several types of members. Finer divisions of a plate require larger 
assembled matrices. Any procedure which would reduce the order of 
these matrices, without increasing the error of the result, would, of 
course, save computer time and space. 

It was anticipated that use of the newer element would allow coarser 
plate divisions with no increase of error. If so, a significant reduc- 
tion in the size of the assembled matrices might be possible in spite of 
the higher order of the element matrices. Hence, the investigation was 
undertaken to see whether the newer element produces a significant im- 


provement. 


PROCEDURE 
The equations for bending energy and work of in-plane forces were 
written in matrix form, making use of the finite element method with 
the sixteen-degree-of-freedom element. A computer program, with matrix 
iteration, was written to generate and solve the matrix equation re- 


sulting from application of conditions of minimum potential energy. 
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Several previously solved cases of plate buckling were run. Results 
were compared with known values and with values obtained using the twelve- 
degree-of-freedom plate element. A few previously unsolved cases were 


investigated. 
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Z. FORMULATION 
INTRODUCTION 
In this section, the formulation of the problem is introduced by 

a discussion of the energy method. The plate element, displacement 
vector, displacement polynomial, and transformation matrix are explained 
as preliminaries. With this background established, the matrix equa- 
tions for bending energy and work of in-plane forces are derived for a 
single plate element. Extension of the method to multiple-element 
plates and the method of handling boundary conditions are discussed next. 
Then, the method of finding a matrix equation for minimum potential energy 
is shown. Finally, the solution of the resulting characteristic value 


problem is discussed. 


THE ENERGY METHOD 


If a plate loaded by external, in-plane forces is bent a small 
amount, the external forces move because the plate edges move. Thus, 
the forces do work. At the same time, to bend the plate, the strain 
energy of bending must be supplied. If the work of external forces 
exceeds the bending energy, the plate buckles. These two energy terms 
are equal in magnitude at the critical load. Or, the total potential 


energy is minimum at the critical load. 


PRELIMINARY DEFINITIONS 


Plate Element and Displacement Vector 





The plate element used is shown above. The nodal order used is 
indicated at the corners of the element. The generalized displacement 


vector used for a single node is* 
Or. = Cw Way Wx Wixy J) (1) 


The displacement vector used for a plate element is 


* 

An underline indicates a matrix. Capital letters are used for 
Square or rectangular matrices. Small letters are used for column 
vectors. A superscript T is used to indicate a transpose. 
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So = C Gyl0,0) ! dplp,0) ! 63(0,q) | S4(pq)I™ (2) 


Displacement Polynomial and Transformation Matrix 


Transverse deflection was represented by 


7 ees 
WED 2 ij x y (3) 


pal j= 


By using this polynomial and derivatives of it, evaluated at the plate 


element nodes, the generalized deflection of the nodes may be represented 


by 
do =Aa (4) 


where A is a 16 x 16 transformation matrix, and a is a vector of the poly- 


nomial coefficients. The polynomial coefficients are then given by 


ag =A'ds (5) 


BENDING ENERGY 


The integral expression for the strain energy of bending “ee 


Wi = 5 EF (Mywyxy + My Wiyy —2MyyWixy dx dy (6) 
This expression may be re-stated in matrix form as 


U = -$Sfml"z dx dy (7) 


il, 


where m and z are vectors defined by 
m=DLMx My -2MyyJ™ (8) 


Z = CWiex Wyy Wry a (9) 
Bending moments are related to curvatures and material properties yo 

My = -D (Woxx +4 Wyyy? 

My =~ D(Woyy + VY Woex) 

My = OU-V) Wy xy 


(10) 


Here, D is the plate flexural rigidity, and v0 is Poisson's ratio. 


Equations (10) may also be written as 


Mx | : v 0 W, xx 
My = -D 4) 4 0 Wryy (11) 
—2 Mxy 0 0 2(1-) Wry 


Defining the matrix of functions of ) to be P allows equation (11) 


to be given in the more compact form 
m=-DPz (12) 


Combining equation (7) and equation (12) gives the new equation for 


bending energy 
U = ZD§f(Pz"Zz dx ady (13) 


Since the vector z is made up of derivatives of the displacement poly- 


nomial, it may be related to the polynomial coefficients by 
Z= Ga (14) 
where G is a 3 x 16 matrix of functions of x and y. Using equation (14) 


to substitute for z in equation (13) gives 


DIS(PGal Ga dx dy (15) 
16 


IN) i 


= 


Substituting for the vector a using equation (5), extracting terms not 


containing x and y from the integrand, and rearranging gives 


U= OS wa P Gdx dy) Ad, (16) 
The portion 
K, = DIAS @ PG dxdy)A" (17) 


is the stiffness matrix for the plate element. Finally, the bending 


energy of the plate element, in terms of the element stiffness matrix, 
is 


U = $ So Kodo (18) 


WORK OF IN-PLANE FORCES 
The integral expression for the additional work done by in-plane 


forces because of lateral deflection iow 
= -} : 2 +2N )dxd 19 
W = -3 Sf (Naw + Ny Wy +2 Nay Wx Wy) dx dy (19) 
The matrix form of this equation is 


W = -3SfrTN* rr dxdy (20) 


Here, r and N* are defined by 


W, 
CoS | (21) 
r Wx 


oe hi Ny | 


aad Nxy Nx —— 





ome 
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The matrix of load parameters may be expressed as 
N*= nN (23) 


where rm is a scalar factor, and N is a matrix of constants which give 
the relative magnitudes of the types of loads. The vector r, which con- 
sists of derivatives of the displacement polynomial, may be related to 


the polynomial coefficients by 
r= Hd (24) 
H is a 2 x 16 matrix of functions of x and y. Equations (20), (23), 


and (24) lead to 


W =- an SJ(Hal"NHa dxdy (25) 


Continuing in the same manner as was used for bending energy formulation 


gives 
W = -5 nde (A*)" (ffHTNH dx dy) Ao (26) 


The portion 
Bo =(A*) SSH NH dx dy) A” (27) 
is the stability coefficient matrix for the plate element.* 
The work of in-plane forces, in terms of the stability coefficient matrix, 
is 


W= -5 150 Bodo (28) 


*The stability coefficient matrix is sometimes referred to as the 
geometric stiffness matrix. 
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EXTENSION TO MULTIPLE ELEMENTS 

So far, the stiffness matrix and stability coefficient matrix for 
a single plate element have been derived. Poor accuracy results when a 
plate is represented by a single element. Rather, the plate must be 
represented by an assembly of plate elements. Although the plate need 
not be divided into identical elements, identical elements were used 
throughout this investigation. 

The same method of extension to multiple elements for the stiffness 
matrix applies to the stability coefficient matrix. Only a brief dis- 
cussion will be given. A more complete description is contained in 
references 7 and 8. 

Let £. be a generalized force vector representing externally ap- 
plied forces and moments at nodes of a plate element. Then, £ and the 


corresponding nodal displacement vector are related by 


fo = Ko do (29) 
where Ko has an order equal to the number of degrees of freedom of the 
plate element. Similarly, when several elements are assembled to make 
up a plate, 

Poe Wwe (30) 
where K is the assembled stiffness matrix, of higher order than K be- 


cause there are more nodes. Elements ofuk are sums of elements taken 


from the plate element stiffness matrices. 


BOUNDARY CONDITIONS 
One of the appealing characteristics of the finite element method is 
the simplicity of handling boundary conditions. For example, no action 
whatever is required for a free edge. For other edge conditions, two ap- 
proaches are available. One method is to eliminate appropriate rows and 


columns of the stiffness and stability matrices, reducing the order of the 
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matrices. For ease of computer programming, a second method was used. 
Zeros were placed in appropriate rows and columns of — matrices. Then 
each zero diagonal element of the stiffness matrix was made unity, be- 
cause the inverse of the matrix was needed later. Table I shows, for 
edge nodes, the nodal displacement components which are zero for the 
prescribed edge condition. It should be noted that the constraints given 
in Table I are the only edge constraints. There is no constraint on in- 
plane deflection for any edge condition; i.e., the plate is "supported" 


or "clamped", but not "held". 


MINIMUM POTENTIAL ENERGY 
Let K and B now represent the assembled stiffness and stability 
matrices, altered by edge conditions. The total potential energy of 
the plate is the sum of bending energy and work of in-plane forces. 


From equations (18) and (28), 


V=U+w=5d' KS - 5nd" Bd (31) 


The critical condition for stability is an equilibrium configuration 
with a set of displacements which makes the total energy a minimum. This 


requires that 
Kd -n Bd =O (32) 


CHARACTERISTIC VALUE PROBLEM 
With some alteration of form, equation (32) expresses a well-known 


characteristic value Seonten 


(DK") 5 =n Bd (33) 
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TABLE I 


BOUNDARY CONDITIONS 


a I 
Edge Condition W Woy W, Wi, 


Simply supported edge 
parallel to x axis 0 0 


Simply supported edge 
parallel to y axis 0 0 


Clamped edge 0 0 0 0 


Line of symmetry 
parallel to x axis 0 0 


Line of symmetry 
parallel to y axis 0 0 


Line of anti-symmetry 
parallel to x axis 0 0 


Line of anti-symmetry 
parallel to y axis 0 0 
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(D/n) K*d = BO | (34) 
AK*S = BS (35)% 


Here, A is an eigenvalue parameter which may have as many values as 
the order of the matrices. Since only the lowest load is desired, 
however, only the largest eigenvalue needs to be found. 


Most plate buckling coefficients are expressed in the cow 


k = (Neiy b7)/ (47? D) (36) 
where b is the plate width in the y direction. If the largest value 
of A is known, the buckling coefficient may be given in the same 


form by 
k =(n b°)/ (PD) = b*/ (at?) (37) 


SOLUTION OF CHARACTERISTIC VALUE PROBLEM 
Solving equation (35) for the largest eigenvalue is the only 
problem remaining. Matrix iteration was the method used, although 
there are other methods. 
The first procedure tried was as follows: 
1. Assume a starting trial vector, 6 . 


2. Calculate an improved trial vector, 6‘, using 


6'= K™ BS (38) 


3. Use the Rayleigh quotient to find an estimate for A _ , using 


A= (8) BS')/( CY KS") (39) 


4. Compute k, using 


*This equation is often given in the form (K*) 1 B-ADd=0. 
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ois / (A 1) (40) 
5. Use the improved trial vector as the next starting trial 

vector, and iterate until successive values of k agree 

within a prescribed limit. 
The procedure did not give good results for shear loads. For pure shear, 
corresponding to each positive eigenvalue, there is a negative eigen- 
value of equal magnitude. The Rayleigh quotient does not give a good 
estimate for A. 


A new procedure was used which proved to be satisfactory. A new 


quotient, 
* =((S'KS') ASKS) (41) 
was used instead of the Rayleigh quotient. Otherwise, the procedure 


was the same as before. 
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3, RESULTS 
INTRODUCTION 
General 
In this section, results of three kinds are given. First, compari- 
sons are made with known solutions to show the accuracy of the method. 
Second, comparisons are made with results given by Kapur and Hartz to 
show the effect of the sixteen-degree-of-freedom element. Third, some 
previously unsolved cases of buckling are given to show the versatility 
of the finite element method. 
Information for Reading Tables of Results 
Wherever appearing, the known values of the buckling coefficient 
are from the methods given by Timoshenko and Gere, but are not necessar- 
ily the values stated by even Some values were computed to a greater 


precision. In all cases, the buckling coefficient is given in the form 
K = (Nerit bY) / (tr? D) (42) 


ACCURACY OF THE RESULTS 
General 
Table II shows the effect of grid size upon the error of the buck- 
ling coefficient for several cases. Table III shows the effect of as- 
pect ratio upon error for a single case. In many cases, the error is 
less than one per cent for rather coarse grid sizes. 


Factors Affecting Earor 


General observations concerning the results are as follows: 


ies Error is greater for coarser grids, as expected. 
25 Error is greater for pure shear loading. 

on Error is greater for clamped edges. 

4, Error is greater for narrower or longer plates. 
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TABLE II 


COMPARISON OF BUCKLING COEFFICIENTS WITH KNOWN VALUES FOR 
SEVERAL CASES, USING SEVERAL GRID SIZES 


Uniaxial Compression 
Aspect Ratio = l 
Poisson's Ratio = 0 


2S k(known) = 4.00000 
SS SS Grid Size k %, Difference 
2 ee 4.01575 +0.4 
a 4.00324 +0.08 
SS a 400104 40.03 
6 x 6 4.00021 +0.005 
8 x 8 4.00007 +0.002 


Isotropic Compression 
Aspect Ratio = l 

2 Poisson's Ratio = 0 
k(known) = 2.00000 


se SS Grid Size k % Difference 
SS 2 az 2.00791 +0.4 
4x4 2.00052 +0. 03 

8 x 8 2.00003 +0.002 


Isotropic Compression 

Aspect Ratio = l 

Poisson's Ratio = 0 

5.30< k(known) < 5.33--lower 
bound used for comparison. 





Grid Size k % Difference 
Zz x 2 5.36223 +1.2 
S23 5.36599 ie 
a eg 5.32713 +0.5 
6 x 6 5.30900 3 Os 
om 6 5.30544 +0.1 
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TABLE II (continued) 











Pure Shear 
Aspect Ratio = l 
—~-S$$ ——— Poisson's Ratio = 0 
| k(known) = 9.34 
SS Grid Size k %, Difference 
5S | Ex Z 10.01581 oye 
3°x" 9.57700 +2.5 
4x4 9.41822 +0 .8 
Uniaxial Compression 
Aspect Ratio = 1 
CL Poisson's Ratio = 0.25 
k(known) = 1.70 
$s oS Grid Size k % Difference 
ae 1.71399 +0.8 
FR 4x4 1.69942 * 
8x4 1.69895 * 


*Zero for precision of known value. 


Pure Shear 
Aspect Ratio = 1 
See 6) |) Poisson's Ratio = 0 
| | k(known) = 14.71 7 
CL cL Grid Size k % Difference 
c 2 x 2 23.26396 +58 .2 
L 3x3 1604630 +9.1 
4x4 15.04271 +2 .3 





9S SS 


FR 


oo 

CL CL 
99 
CL 

35 98 
CL 


TABLE II (Continued) 


Uniaxial Compression 
Aspect Ratio = 1 


Poisson's Ratio = 0.25 


k (known) = 1.43418 


Grid Size k 
en ae 4 1.44324 
ay & 1.43483 
8x4 1.43435 
Pure Shear 


Aspect Ratio = l 
Poisson's Ratio = 0 
k(known) = 12.28 


Grid Size k 
25X22 17.13879 
4m 4 12.82566 


Uniaxial Compression 
Aspect Ratio = 1 
Poisson's Ratio = 0 
k(known) = 6.74 


Grid Size k 
2&2 6.82334 
Axx 4 6.78242 
8 x 8 6.74583 


Uniaxial Compression 
Aspect Ratio = l 
Poisson's Ratio = 0 
k(known) = 7.69 


Grid Size k 
Jace? 8.79292 
4x4 7.74665 
8 x 8 7.69541 


2/ 


% Difference 


+0 .6 
+0 .05 
+0.01 


%, Difference 


+39.6 
+4.4 


% Difference 


ol  /4 
+0 .6 
+0 .09 


% Difference 


+14 .3 
Hohe? 
+0.07 


TABLE IIT 


COMPARISON OF BUCKLING COEFFICIENTS WITH KNOWN VALUES 
FOR ONE CASE, USING SEVERAL ASPECT RATIOS 








Ss na Uniaxial Compression 
<j a Grid Size, 4 x 4 
= Poisson's Ratio = 0 

b oh Aspect ratios ror byck}ing.mede 
\ transition: 2%, 3%, 42, etc. 

| oD 

: a/b k k (known) % Error No. Half Waves 

O22 27.05296 27.04000 +0.048 I! 
0.4 8.41330 8.41000 +0 .039 1 
0.6 5.13940 5.13778 +0 .032 a 
0.8 4.20364 4.20250 +0.027 1 
1.0 4.00104 4.00000 +0.026 il 
dvs 4.13556 4.13434 +0 .030 1 
1.4 4.47150 4.47022 mo. O2¢ iL 
ZO 4.00770 4.00000 +0.193 2 
30 4.02924 4.00000 +0. 7 3 


When the buckled shapes of the plates are considered, all of these 
observations may be further generalized to a single one. The accuracy 
of representing the buckled surface is less for greater curvature gradi- 
ents. To keep error low, finer grids must be used for more complex 
buckled shapes. 

Sources of Error 

Discretization Error Most numerical methods replace continuous 
equations with discrete approximations. The finite element method, 
however, uses continuous equations but a discrete approximation of the 
continuous structure being analyzed. The finite element method has in- 
herent discretization error. Since the finite element approximation 
of a real plate is stiffer than the real plate, the discretization 
error of the buckling coefficient is never negative. In addition, the 
discretization error decreases montonically as the plate is further 
subdivided, provided that previously existing nodes are also nodes in 
the finer weiae> For errors shown in Tables II, III, and IV, most 
error is believed to be discretization error. 

Manipulation Error Matrix iteration produced monotonically con- 
verging values of the buckling coefficient. Further, as a consequence 
of the matrix iteration method used, convergence was always from an 
initially higher value. Hence, matrix iteration error of the buckling 
coefficient is also positive. A convergence criterion of .001 per cent 
difference between successive values of the buckling coefficient was 
used. Consequently, the error of the results cannot be much less than 
.OO1 per cent. Round-off errors are believed to be insignificant in 
comparison with other errors, since double precision programming was 


used. 
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Negative Error 


Because the two principal sources of error--discretization and 
convergence--both produce positive error of the buckling coefficient, 
the error cannot be negative. Although some errors initially appeared 
to be negative, in each case encountered the cause was found to be lack 
of precision of the "known" value of the buckling coefficient. 

Accuracy of Buckling Mode Vectors 

Although the buckled shape of the plate is of secondary interest, 
the buckling mode vector was available as a consequence of the method used. 
The vector provided another check on the accuracy of the method. For 
example, the nodal displacements were compared with values computed from 
the known buckled shape for the uniaxially compressed, simply supported 
plate of aspect ratio one for a 4 x 4 grid. The displacements were in- 


ternally consistent to a minimum of three significant figures. 


EFFECT OF SIXTEEN-DEGREE-OF -FREEDOM ELEMENT 

Table IV shows comparisons of results with those given by Kapur and 
Hartz for two cawese An improvement in accuracy is apparent. More 
important, for the same order of magnitude of the error, fewer elements 
were required. 

To show what this improvement means, consider the case of the simply 
supported plate compressed in one direction. Kapur and Hartz report an 
error of -0.58 per cent when using a 12 x 12 grid. Assuming that no 
symmetry and no compaction of matrices is used, and not considering re- 
duction of the matrix order when boundary conditions are applied, their 
formulation requires a 507 x 507 assembled stiffness matrix to represent 
three degrees of freedom at each of 169 nodes. A stability coefficient 


matrix of the same order is also required, of course. Using one more 
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TABLE IV 


COMPARISON OF BUCKLING COEFFICIENTS WITH VALUES OBTAINED 
USING TWELVE -DEGREE-OF-FREEDOM ELEMENT@ 


3S Uniaxial Compression 
Aspect Ratio = l 
— Poisson's Ratio = 0 


oS SS k(known) = 4.00000 
SS 12 Degrees of Freedom 16 Degrees of Freedom 
Grid Size k % Error k % Error 
2 KkKKK kkkKK 4.01575 +0 .39 
3 °x 3 3.645 -8 .88 4.00324 +0 .08 
4x4 3/8 -5.75 4.00104 +0 .03 
6 x 6 3.887 -2.80 4.00021 +0.005 
8 x 8 3.933 -1.68 4.00007 +0.002 
10 x 10 3.96 -1.00 RKKKKKK KEKKEK 
12 x 12 3.977 -0.58 KEKKKKE KRKKKEK 
CL Isotropic Compression 
Aspect Ratio = l 
Poisson's Ratio = 0 " 
CL CL 5.30< k(known )< 5.33 
12 Degrees of Freedom 16 Degrees of Freedom 
CL 
Grid Size k % Disget . k Ze Diet . 
2a 2 4.901 -7.53 5.36223 +1.17 
3x3 KKK KKKKK 5.36599 +1.24 
4Gx& 4.975 -6.13 5.32713 +0.51 
6 x 6 5.078 -4.19 5.30900 +0.17 
8 x 8 5.160 -2.64 5.30544 +0.10 
10 x 10 5.2 -1.57 RKKKKKK KEKKK 


“pata from Kapur and Hartz, reference 4. 


Othe lower bound was used for computing per cent differences. 
Since Kapur and Hartz had used an average of the upper and lower bounds, 
the per cent differences given by them were re-computed. 
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degree of freedom per node, the error is 0.39 per cent when using a 

2x 2 grid. For the same assumptions, the stiffness matrix size is 36 x 
36 for four degrees of freedom at each of nine nodes. In this example, 
the superior accuracy of the newer plate element and the consequent 


savings in computer time and space are unequivocally apparent. 


PREVIOUSLY UNSOLVED PROBLEMS 

Table V shows the results of investigation of two previously un- 
solved plate buckling problems. 

An edge which has more than one edge condition along its length 
is difficult to analyze by other means. Although the computer program 
used for other cases had to be altered very slightly to handle such an 
edge, there was no particular difficulty in using the finite element 


method. 
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TABLE V 


PLATE BUCKLING PROBLEMS NOT PREVIOUSLY INVESTIGATED 


FR 





Uniaxial compression, aspect ratio = 2, Poisson's ratio = 0.3 


Nodes marked X are either all simply supported or all clamped. 


—_elceElUOEeTetlClCC rE lUlUlCU ST ll eee eee eee eee eee eee less eee tl eee eee ll ees ess ees ess sss ese esse eee eee eee eee eee eee eee eel eel eee eee 


HALF FREE, HALF SIMPLY SUPPORTED 
k = 0.179 


Relative Transverse Displacements of Nodes 


Node W Node WwW Node W Node Ww Node W 
1 -.3562 2 -.1207 5 0 4 0 5 0 
6 -.3662 7 -.1354 8 -.0139 9 +.0019 10 0 
ll -.3699 12 -.1408 13 -.0182 L4rr=F90027 15 0 
16 -. 3662 17 -.1354 18 -.0139 19 +.0019 20 0 
21 -.3562 22 -.120/7 26 0 24 0 25 0 


HALF FREE, HALF CLAMPED 
k = 0.208 


Relative Transverse Displacements of Nodes 


Node W Node WwW Node W Node WwW Node Ww 
1 -.3174 jd -.0956 3 0) 4 0 5 0 
6 -.3257 7 -.1073 8 -.0034 9 +.0013 10 0 
11 -.328/7 £2 -.1111 13 -.0053 14 +.0026 15 0 
16 =, 329/ 17 -.1073 18 -.0034 19 +,.0013 20 0 
Pad -.3174 22 -.0956 23 0 24 0 25 0 








4. THE COMPUTER PROGRAM 
GENERAL 
A sample computer run is given in Appendix I. The program listing 
is largely self-explanatory. The program is in FORTRAN IV with double 
precision. The IBM 360 computer at the Naval Postgraduate School was 
used. This computer has a core memory capacity of 512,000 bytes. The 


HASP compiler used required about 55,000 bytes. 


CAUTIONS 

Program Efficiency 

The program was written as an aid in producing this thesis only, 
by a relatively inexperienced programmer. Paradoxically, one purpose 
of the thesis is to show how computer time and space can be saved, yet 
both were often sacrificed for ease of programming. 
Convergence in Matrix Iteration 

Initially, a trial vector with a high content of some mode shapes 
was used to start the matrix iteration process. This sometimes caused 
convergence to an incorrect mode. So that the trial vector would not 
be prejudiced toward particular modes, a new trial vector consisting of 
elements taken from a table of random numbers was used. In all cases 
run using the new trial vector, where the proper buckling mode was 
known, convergence was to the correct mode. 

Convergence was very slow for cases with the two largest eigen- 
values nearly equal. Although methods of accelerating convergence were 
tried, no generally acceptable method was found. The program is not 


recommended for such cases. 
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5. SUMMARY 
CONCLUSIONS 
For plate buckling problems using the finite element method, 
the sixteen-degree-of-freedom element is far more accurate than the 
twelve-degree-of-freedom element. Not only is the accuracy superior, 
but it is sufficiently greater so that sizes of assembled stiffness 
and stability coefficient matrices may be greatly reduced. Far fewer 
plate elements are needed for equivalent accuracy. Computer time and 


Space requirements are significantly reduced. 


PROBLEMS NEEDING FURTHER STUDY 
The following problems were beyond the scope of this investigation: 
ibe An attempt should be made to find a general rule for sub- 
dividing the plate for minimum error. Elements of un- 
equal sizes, or more divisions in one direction might be 


used, for example. 


Vig A systematic investigation of many of the remaining unsolved 
cases of plate buckling should be made. 


3, Further study of methods of determining the eigenvalue 
should be made in order to find a more rapid method. 


4. A more general and more efficient computer program 
should be devised for industrial use. 
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APPENDIX I 


SAMPLE COMPUTER RUN 
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1 
ESTABLISH TRIAL VECTOR OF FLEMENTS RETWEEN -1 AND ¢1 TAKEN FROM 
TABLE OF RANONM OIGITS 
II=1 
JJ=10 

497 DO 739 I=1,10 
READ (5,740) (WOJ),J=I1,JJ) 
1I=11410 

= Sd 

146 HURMAT Viogs.2) 

741 FORMAT (1HO, LOHTRIAL VECTOR, FIRST,14,1X,17HELEMENTS OF ARRAY) 
JJ=10 
DO 498 I=1,10 
WRITE (69742) (W(J) eJ=II 2 JU) 

Il=1f*1lo 
JJ=JJ+10 

498 CONTINUE 

742 FORMAT (10F6.2) 

MULTIPLY W TRANSPOSE * S * W = WTSW. 

700 wISwW=0. 
00° 701 I=1, INDEX 
00 701 J=1, INDEX 

701 WISW=WISwew(1) *STIFF( I,J) *Ww(J) 

NORMALIZE We 


702 


611 
613 
612 


qh) 


NAIN 
C00 
aw 


706 


606 


708 


715 


714 


DO 702 IT=1,I1NDEX 
WOT) =WCT)/DSORTCWT SW) 
IF (MM.ZEQ.1) GO TO 713 


AFTER FINAL ITERATION,WRITE TRANSVERSE NISPLACEMENTS IN MAP FORMAT 


IF (C3.EQ.0.) GO TO 721 

DO 716 I=1,fNDEX 

[IF (MM.EQ.2) WITD=WSAVECT) #W(T) 

IF (MM.EQ.3) WC )=20*#WSAVE(I)-W(T) 
CONTINUE 

WRITE (6,608) 

FORMAT (1H1,22X»24HTRANSVERSE DISPLACEMENTS) 
WRITE (6,609) 

FORMAT (1HO,24X,20HY INCREASES TO RIGHT/) 
WRITE (67610) 

FORMAT (21X,26HX INCREASES TOWARDS BOTTOM/) 
L=4*N-3 

DO 612 J=1,M 

WRITE (69611) (WOUI) pI=Kyl > 4) 

FORMAT (8(011.4,5X) 3 

WRITE (6,613) 

FORMAT (1H S/S/SSISIS) 

K=L+4 

L=L+4eNn 

IF (C3.£0.0.) GO TO 712 

MM=MMel 

IF (MM.EQ.3) GO TO 722 

IF (MM.E0.4) GO TO 712 


MULTIPLY SINV * 8 * WI NORMALIZED) = WNEXT 
Osos T=1, INDEX 


MULTIPLY WNEXT TRANSPOSE * S * WNEXT = EIGEN 


EIGEN=0. 
DO 706 I=1,INDEX 
I oe 


DQ 706 J=l, 

EIGEN=E IGEN¢WNEX TCL) *STIFE( I,J) *WNEXT(Q) 
FIND BUCKLING COEFFICIENT. 

AK=]./(DSQRT (EIGEN) *(3.14159265*AR) ¥#2) 
wRItE {2 2665 L AK 

FORMAT (1HO,9HITERATION,1394H, K=,F9.5) 
TEST FOR 20 ITERATIONS. 

IF (L.EQ.20) GO TO 710 

TEST FOR CONVERGENCE. 


IF (L-EQ.1) GO TO 708 
Be ee OABSTSAVE) 7048S TAROT <0 60 TN 711 
L=L+¢l 


SAVE W (NORMALIZED) 


DO 715 I=l,INDEX 
WSAVE(I)=WIT) 


USE RESULTING VECTOR AS NEXT TRIAL VECTOR. 
00 709 T=1,INNEX 
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TO9 WCT)=WNEXTOT) 
GA TN 700 

T10 WRITE (6,522) 

522 FORMAT ({HO, 33120 TTEFRATIONS WITHOUT CONVERGENCE) 
GC TO 712 

TLL WRITF (6,525) 

Bese. cell (LHO,48HSNICKL ING PARAMFTFR HAS CONVERGED TN .NOONL PFRCFMT) 
MM=? 
GO ™ 722 

T1? CONT INUF 
FANN) 

ASPECT RATIO= 1.00000 POISSON RATIN=0.0 NY. ELEMENTS IN X=4 NO. FLEMFNTS IN Y=4 


AGUNDAPY CONDITIONS, L=FREE, 2=PINNEID, 3=CLAMPED, 4=SYM., S=ANTI-SY™. 
Y=O,EVGF CONDITION 3 Y=8,FNGE CONOITION 4 X=O0,FOGE CONDITION 2? X=A4,FNGE CANNTTIAN © 
X<FLATTIVE LOAN ST7ES NX= 1.00 NY= 0.0 NXY= 0.9 


TRIAL VECTOR, FIRST : NTS AR 
0.0 hoS9 +f 485 0129 Seal 8.55 aOR ay 0.68 -0.54 0.64% 
—0.09 -Oste=0F65 O.im O25" Off) -0O.17 -0.594 .72 99.75 
-0.09 -Q.11 0.60 -0.79 On 7) -0.C9 O. 63 0-18 -0.17 G.58 
Occ 0.14 0.91 0.14 O12 Oak OeS6 -0.93 OP 8) -0.43 
O.1s -0.47 Oem -O.75 -0.55 O.AD Oeste? 0.98 -0.25 0.04 
0.37% 0.42 -0.98 0.97 -0.60 -0.04 -0.18 0.30 0.15 -0.4! 
0-31 -0.65 -0.57 -0.75 -0.76 0.97 0.21 9.98 -0.21 0.2) 
O74 -C. 36 0.66 -0.68 -0.19 -0.4°9 -0.76 0.55 ~Ge 33 —@.. Aid 
0.53 0.34 0.99 0.95 0.38 -(). 39 0.41 -0.95 0.98 o.tS 
et | -~OR 74 0.34 0.99 Oe -0.31 0.94 -0.05 -0.459 -0.23 


ITERATION Ly K=104.5164617 
ITERATION 2,5 K= 26.5343? 
[TFRATICN 3, K= 13.0965) 
ITERATION 4y K= 8.81542 
ITFRATION 5S, K= 7.89034 
[TEFRATIGN 6, K= 7.72830 
ITFRATICN V7, K= 7.270094 
ITERATION 8, K= 7.99634 
PYERATIC™N S, K= 4.69555 
ITERATION 10, K= 7.09543 
ITFRATIGN Ll, K= 7.69541 
SUCKLING PARAMETER HAS CONVERGEN TO .001 PERCENT 


4535 


TRANSVERSE OF SPLACEMENTS 


Y INCREASES TO RIGHT 


X INCREASES TOWARDS BOTTOM 


0.0 


-0.1212D0-01 


-0.17150-01 


-0.1214D-01 


0.0 


-0.3320D-01 


-0.46970-01 


-0.33230-01 
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-0.5004D-01 


-0.70810-01 


-0.5009D-01 


=) ...0320-00 


-0.79690-01 


-0.5538D0-01 
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